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Introduction 


There  are  always  interesting  problems  in  engineering  and  science  that  are  too  big  to  fit  on  available 
computers.  Regardless  of  whether  an  analyst  is  using  a  personal  microcomputer  or  a  large 
supercomputer,  the  issue  of  available  memory  is  an  important  one  that  often  limits  the  type  and 
scope  of  analysis  that  is  to  be  performed.  While  increases  in  density  of  computer  memory  and 
efficient  vinual  memory  management  schemes  have  helped  to  alleviate  the  "supply"  side  of 
memory  requirements,  reduced  coordinate  methods  (commonly  termed  "modal"  methods)  arc  an 
important  means  to  help  reduce  the  "demand"  for  more  memory.  These  reduction  schemes  can  be 
thought  of  as  techniques  to  "make  large  problems  small",  while  retaining  as  many  of  the  important 
aspects  of  the  solution's  behavior  as  possible.  Because  a  large  problem  is  noi  small,  reduced 
coordinate  algorithms  have  the  aura  of  "getting  something  for  nothing",  and  so  are  an  appealing 
area  for  study  in  computation.  This  report  attempts  to  document  some  sample  benefits  and  risks  of 
using  these  modal  analysis  procedures  in  Civil  Engineering  applications. 

The  application  of  modal  methods  to  large  problems  in  Civil  Engineering  has  a  rich  history.  The 
classical  technique  of  "normal-mode  integration"  from  earthquake  engineering  is  an  example  of  a 
modal  superposition  method  that  has  been  used  with  great  success  on  a  wide  variety  of  linear 
problems.  Extensions  of  this  approach  to  nonlinear  problems  are  much  more  difficult  firom  the 
dual  standpoints  of  theory  and  application.  Strictly  speaking,  the  principle  of  superposition  is  not 
valid  in  a  nonlinear  environment,  but  appeals  to  local  linearization  of  the  underlying  problem  can  be 
used  to  extend  the  mode  superposition  principle  to  this  case.  In  a  linear  problem,  one  set  of 
representative  modes  can  be  calculated  at  the  outset,  and  used  throughout  the  entire  temporal 
analysis.  Unfortunately,  the  effect  of  nonlinearities  in  a  problem  may  be  to  turn  an  initially 
satisfactory  set  of  modes  into  a  poor  estimate  of  the  actual  response  of  a  system  once  substantial 
nonlinearities  have  developed.  Development  of  robust  schemes  to  update  or  augment  a  set  of 
modes  to  account  for  solution  nonlinearitics  is  presently  an  important  unanswered  research  issue 
that  must  be  solved  for  the  widespread  use  of  nonlinear  modal  analysis  to  occur. 

The  theory  and  practice  of  modal  analysis  techniques  are  presented  in  Chapter  2  in  the  mathematical 
framework  of  projection  methods.  The  example  applications  of  Chapter  3  show  the  results  of 
application  of  reduced  coordinate  schemes  to  a  few  representative  problems  in  Civil  Engineering. 
The  final  chapter  reiterates  the  results  of  the  applications,  and  suggests  avenue  for  future  research 
in  this  field. 
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Review  of  Previous  Work 

Modal  projection  methods  have  been  used  in  a  wide  variety  of  solutions  to  problems  in 
engineering,  especially  those  involving  structural  dynamics  and  other  time-dependent  phenomena. 
The  underlying  theory  for  such  projection  schemes  is  developed  in  a  previous  report  [1],  as  is  a 
sample  bibliography  for  this  topic.  Therefore,  only  the  most  important  features  of  these  methods 
will  be  presented  here. 

The  Projection  Theorem 

The  mathematical  setting  for  projection  solutions  is  a  Hilbert  Space,  which  is  defined  as  a  normed 
linear  vector  space,  with  the  norm  (length  function  for  vectors)  induced  by  an  inner  product.  The 
most  important  interpretation  of  a  Hilbert  Space  is  that  these  abstract  mathematical  entities  are  the 
topological  generalization  of  ordinary  three-dimensional  space.  This  means  that  such  concepts  as 
angle,  perpendicularity,  and  completeness  are  preserved,  even  in  an  infinite-dimensional  setting. 

The  fundamental  properties  of  a  Hilbert  Space  are  defined  in  terms  of  the  associated  inner  product 
of  two  vectors.  The  most  important  result  for  seeking  approximate  solutions  to  operator  equations 
defined  on  a  Hilbert  Space  is  the  Projection  Theorem,  which  can  be  used  to  judge  the  quality  of  an 
approximation  family.  The  geometric  interpretation  of  the  projection  theorem  for  ordinary  three- 
dimensional  space  is  shown  in  Figure  1. 


Figure  1 :  Geometric  Interpretation  of  the  Projection  Theorem 


Modal  Projection  Methods 


Page  3 


Given  a  Hilbert  Space  H  (in  this  case,  R3),  and  a  complete  manifold  (subspace)  M,  the  Projection 
Theorem  states  that  every  vector  u  in  H  can  be  uniquely  decomposed  into  the  sum  u  =  v  +  w, 
where  the  vector  v  is  in  M  and  w  is  perpendicular  to  M.  In  the  context  of  approximation,  the 
manifold  M  represents  a  family  of  approximate  solutions,  and  the  vector  w  =  u  -  v  represents  the 
error  in  approximating  u  by  any  vector  v  in  M.  In  this  case,  the  Projection  Theorem  implies  that 
the  length  of  any  error  vector  w  =  u  -  v  can  be  minimized  by  choosing  w  to  be  perpendicular  to  the 
approximation  u,  so  that  v,  w,  and  u  form  an  abstract  generalization  of  a  right  triangle.  The 
notions  of  length  and  perpendicularity  depend  on  the  underlying  definition  of  the  inner  product  that 
defines  the  Hilbert  Space. 

Two  concrete  examples  of  application  of  the  Projection  Theorem  are: 

1.  Galerkin  Finite  Element  approximations  for  self-adjoint  boundary-value  problems.  In  this 

case,  the  Hilbert  Space  is  the  space  of  all  admissible  solutions,  and  the  inner  product  function 
is  some  type  of  associated  potential  energy  functional  for  the  boundary-value  problem.  The 
approximation  manifold  is  the  span  of  the  Finite  Element  basis  functions,  and  the  error  vector 
is  the  usual  residual  function.  The  familiar  error-minimization  interpretation  of  Galeridn's 
criterion  for  choosing  the  best  approximate  solution  from  the  family  is  the  minimum  error 
vector  length  guaranteed  by  the  Projection  Theorem. 

2.  Reduced  coordinate  approximations  for  the  Finite  Element  equations  resulting  from  Example  1. 

In  this  setting,  the  underlying  Hilbert  Space  is  the  fmite-dimensional  space  of  all  nodal 
solution  vectors,  and  the  approximation  manifold  is  the  range  of  all  matrices  whose  columns 
are  the  modes  used  for  approximation.  The  inner  product  can  be  the  "usual"  dot  product  for 
vectors,  or  a  weighted  product  induced  by  multiplication  by  the  Finite  Element  mass  or 
stiffness  matrices. 

In  this  setting,  it  appears  that  there  is  a  distinct  analogy  between  the  following  two  questions: 

Ql.  How  many  nodes  (or  basis  functions)  are  required  to  obtain  an  acceptable  solution  quality 
for  a  Finite  Element  Model,  and  how  should  these  nodes  be  chosen? 


Q2.  How  many  modes  are  required  to  obtain  an  acceptable  reduced  coordinate  approximation 
for  a  set  of  Finite  Element  equations,  and  how  are  these  modes  to  be  constructed? 
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Other  standard  concerns  involving  Finite  Element  approximation  families  (such  as  the  use  of  mesh 
rezoning  to  update  the  solution  family  after  permanent  deformations  have  occurred)  also  have 
parallel  interpretations  in  the  setting  of  reduced  coordinate  models. 


Eigenvector  Expansion  Methods 

While  the  Projection  Theorem  can  be  used  to  choose  the  best  approximation  from  a  family  of 
approximate  solutions,  the  topic  of  the  optimal  choice  (or  construction)  of  this  family  is  one  not 
addressed  directly  by  the  Projection  Theorem.  Insight  into  this  choice  can  be  obtained  by 
considering  the  spectrum  of  the  operator  that  defines  the  problem  to  be  solved. 

The  traditional  approach  to  modal  approximation  for  self-adjoint  problems  has  been  the  use  of 
exact  eigenvectors  as  a  reduced  coordinate  basis.  The  familiar  "separation  of  variables"  method  for 
finding  solutions  to  partial  differential  equations  in  terms  of  (truncated)  Fourier  Series  is  an 
example  of  the  use  of  this  approach,  as  is  the  "normal  mode  integration"  procedure  commonly  used 
in  dynamic  structural  analyses.  Although  there  are  alternatives  to  exact  eigenvectors  that  arc  easier 
to  construct,  the  use  of  eigenvector  expansions  is  a  useful  theoretical  tool  for  consideration  of  the 
choice  of  modes  for  a  reduced  coordinate  analysis. 


If  the  desired  operator  equation  to  be  solved  via  a  reduced  coordinate  method  can  be  written  in  the 
general  form  Ax  =  b  (where  A  is  a  self-adjoint  linear  operator),  then  a  consideration  of  the  quality 
of  various  eigenvector  projection  schemes  can  be  made  by  studying  the  set  of  eigenvalues  of  A. 
This  set  is  known  as  the  spectrum  of  A,  and  when  A  is  a  finite-dimensional  operator  (i.e.,  an  NxN 
matrix),  the  spectrum  consists  of  N  real  numbers  (which  are  not  necessarily  distinct).  For  most 
stable  structural  systems,  these  real  numbers  are  all  positive,  and  in  this  case  A  is  termed  a  coercive 
operator,  or  a  positive-definite  matrix.  A  picture  of  the  spectrum  of  A  for  Uiis  situation  is  shown  in 
Figure  2,  where  A's  eigenvalues  are  diagrammed  in  the  complex  plane. 
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Figure  2:  The  Spectrum  of  a  Positive-Definite  Operator 
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For  structural  dynamics  applications,  the  eigenvalues  Xi  have  the  physical  interpretation  of  squared 
frequencies,  and  so  the  minimal  eigenvalues  are  associated  with  the  lowest  frequencies  of 
vibration,  and  the  maximal  ones  with  higher  frequency  effects.  When  the  frequencies  of  a 
discretized  system  (such  as  a  Finite  Element  model  of  a  structural  system)  are  compared  with  the 
exact  frequencies  of  the  actual  continuous  system,  it  is  generally  true  that  the  lower  approximate 
frequencies  are  more  accurate  estimates  for  the  true  values  than  the  higher  frequencies.  For 
example,  the  fundamental  frequency  for  a  Finite  Element  model  of  a  structure  is  typically  a  very 
accurate  estimate  of  the  fundamental  for  the  actual  structure,  but  the  agreement  between  higher 
frequencies  may  be  very  poor.  Because  of  this  phenomena,  it  appears  that  the  use  of  eigenvector 
expansions  for  approximation  purposes  should  begin  with  the  lowest  frequencies  in  order  to 
achieve  highest  accuracy  for  a  given  number  of  modes. 

Consideration  of  the  expansion  of  the  solution  to  Ax  =  b  in  terms  of  A’s  eigenvalues  lends  further 
credence  to  this  approximation  strategy.  If  the  eigenvectors  of  A  are  denoted  by  ui,  then  the 
solution  X  and  th**  data  b  can  be  represented  in  terms  of  the  eigenvectors  by: 

N  N 

i  =  1  i  =  I 

With  this  convention,  the  unknown  coefficients  that  define  the  solution  x  can  be  found  by  use  of 
the  defining  relation  for  A's  eigenpairs: 


N  N  N 

Aui  =  XiUi  impliesthat  Ax  =  A  ^iUi=  b  =  ^  ^  PjUj  =  b 

i  =  1  i  =  1  i  =  1 

Identification  of  each  component  in  this  equation  gives  a  formula  for  the  coefficients 

Thus  the  solution  components  are  most  sensitive  to  the  smallest  eigenvalues  (the  lowest 
frequencies),  and  the  observation  that  the  eigenvectors  associated  with  the  minimal  frequencies  will 
(in  general)  form  the  best  approximation  is  verified.  This  is  why  the  normal  mode  integration 
schemes  are  widely  used  for  linear  dynamics  problems:  the  eigenpairs  associated  with  the  lowest 
frequencies  lead  to  more  efficient  modal  projection  schemes. 
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Krylov  Subsoace  Methods 

Construction  of  exact  eigenvectors  involves  finding  an  orthonormal  set  of  vectors  (i.e.,  a  new 
coordinate  system)  that  diagonali2es  the  operator  A,  so  that  the  coupled  operator  equation  set 
resolves  into  uncoupled  scalar  forms  when  referred  to  this  new  coordinate  system.  Other 
"distinguished  orthonormal  sets"  exist  that  allow  the  operator  A  to  take  special  simple  forms,  and 
these  alternative  coordinate  systems  may  be  considerably  easier  to  construct  than  a  set  of  exact 
eigenvectors.  Probably  the  most  important  of  these  orthonormal  sets  is  that  generated  by  the 
Krvlov  Sequence: 

Ay,  A2y,  A3y,  A^y,  . 

The  manifold  spanned  by  vectors  from  this  sequence  is  termed  the  Krylov  Subspacc  K(A,y,m) 
K(A,  y,  m)  =  span  (Ay,  A2y,  A3y, . .  A^y) 

It  can  be  shown  [Reference  10]  that  for  self-adjoint  A,  the  vectors  of  the  Krylov  Sequence  form  a 
three-term  recurrence,  so  that  A  is  tridiagonal  when  referred  to  an  orthonormal  basis  for  K(A,y,m). 
This  fact  forms  the  basis  for  the  Lanczos  Algorithm,  which  constructs  this  basis  in  a 
computationally  efficient  manner. 

Since  a  tridiagonal  form  of  the  operator  A  is  not  appreciably  more  difficult  to  use  than  a  diagonal 
one,  there  is  a  strong  incentive  to  use  the  Krylov  Subspace,  since  construction  of  this  manifold  is 
considerably  simpler  than  construction  of  an  equal  number  of  exact  eigenvectors. 

It  is  easy  to  show  that  successive  vectors  in  the  Krylov  Sequence  tend  towards  the  eigenvectors 
associated  with  the  maximal  eigenvalues  by  considering  the  expansion  of  the  initial  vector  y  in 
terms  of  the  eigenvectors  of  A. 

N 

If  y=  the  Krylov  sequence  Ay,  A2y,  A3y .  Amy  can  be  written  as: 

i  =  1 

N  N  2  N  3  ^  m 

. 

i=l  .=1  i=l  i=l 
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As  m  increases,  the  terms  involving  Xn,  Xn-i,  Xn-2,  ...  dominate,  since  these  eigenvalues  are  the 
largest  in  magnitude.  This  leads  to  the  following  convention  (which  is  not  standard,  but  will  be 
used  in  this  document); 

A  Forward  Krylov  method  for  approximating  the  solution  of  Ax  =  b  is  one  that  uses  the  Krylov 
Subspace  K(A,  y,  m)  associated  with  the  operator  A  as  an  approximate  solution  manifold. 

A  Backward  (or  Inverse)  Krylov  method  is  one  that  uses  the  Krylov  Subspace  K(A-i,  z,  m) 
associated  with  the  inverse  of  A  as  an  approximate  solution  manifold. 


Since  the  eigenvalues  of  A*!  are  the  reciprocals  of  those  for  A,  a  backward  Krylov  scheme  uses 
information  from  the  eigenpairs  at  the  low-frequency  end  of  the  spectrum,  and  a  forward  scheme 
uses  data  from  the  high-frequency  end  of  the  spectrum.  It  is  apparent  from  the  discussion  above 
that  the  backward  schemes  are  typically  more  efficient  in  terms  of  number  of  modes  required  to 
form  an  approximate  solution,  but  the  issue  of  expense  of  forming  each  mode  should  also  be 
considered.  Each  vector  added  to  the  forward  sequence  requires  a  matrix  multiplication  by  A, 
compared  to  the  solution  of  a  set  of  equations  involving  the  coefficient  matrix  A  required  for  each 
vector  in  the  backward  sequence.  Thus  each  vector  in  a  backward  scheme  requires  considerably 
more  work  to  form,  which  reduces  the  efficiency  of  the  method.  In  addition,  if  A  is  singular  (or 
nearly  so),  the  solution  of  the  system  Ayi  =  yi-i  to  obtain  the  next  vector  in  the  sequence  may  not 
be  feasible.  This  would  be  the  case  if  the  structural  system  of  interest  possessed  some  type  of 
instability  (another  source  of  near-singularity  will  be  presented  in  the  next  chapter). 

Based  on  this  elementary  analysis  of  Krylov  Subspace  methods  (i.c.,  "all  other  things  being 
equal"),  the  following  general  characteristics  of  the  two  strategies  can  be  contrasted: 

Forward  schemes  will  require  more  vectors  to  achieve  a  given  accuracy  of  approximate  solution, 
but  each  vector  in  the  spanning  sequence  will  be  easier  to  form.  The  approximate  solution 
manifold  will  be  dominated  by  information  from  the  high-frequency  end  of  the  spectrum, 
which  is  less  desirable  from  the  standpoint  of  accuracy,  but  zero  eigenvalues  (corresponding 
to  structural  instability)  should  cause  little  problem  for  this  family  of  methods. 

Backward  schemes  will  require  fewer  vectors  for  given  accuracy,  but  each  vector  will  be  more 
difficult  to  form.  The  solution  manifold  will  be  dominated  by  low-frequency  data  (which  is 
desirable  for  accuracy),  but  the  effects  of  zero  eigenvalues  must  be  alleviated. 


Modal  Projection  Methods 


Page  8 


Examples  of  forward  Krylov  projection  schemes  include  the  Conjugate-Gradient  Method  (CGM) 
[3,8],  and  the  Lanezos  Algorithm  [1,2,9,12,13,14,15].  In  a  nonlinear  time-dependent  problem, 
these  schemes  involve  reconstructing  an  appropriate  Krylov  Subspace  at  each  solution  step  (i.e.  for 
each  equation  of  the  form  A(x,t)  x(t)  =  b(x,t)  that  is  to  be  solved).  This  reconstruction  is  feasible 
only  because  each  vector  in  the  sequence  is  easy  to  form  via  a  matrix  multiplication,  and  has  the 
advantage  that  the  approximation  subspace  "tracks"  the  evolution  of  the  operator  A(x,t)  through  the 
effects  of  both  time  and  sclution  nonlineariiies.  The  biggest  disadvantage  of  these  schemes  is  that 
when  A  is  ill-conditioned  (i.e.,  the  range  of  scales  of  the  frequencies  is  large),  these  methods  may 
converge  very  slowly,  or  not  at  all.  This  difficulty  can  often  be  ameliorated  by  compressing  the 
spectrum  by  premultiplication  of  Ax  =  b  by  an  operator  chosen  to  improve  the  condition  of  A. 

Such  a  matrix  is  termed  a  preconditioner,  and  the  efficient  construction  and  application  of 
preconditioners  forms  the  basis  for  virtually  ail  of  the  practical  Conjugate  Gradient  schemes  [8,1 1]. 

An  example  of  a  backward  Krylov  method  is  the  Inverse  Lanezos  Algorithm,  which  is  used  in  the 
next  chapter  to  find  approximate  solutions  of  some  representative  engineering  problems.  This 
approach  is  the  generalization  of  "normal  mode  integration"  techniques  to  Krylov  approximations 
for  nonlinear  problems.  The  details  of  this  method  have  been  elucidated  in  a  previous  report  [10]. 


Example 

Applications 
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Example  Applications 

This  section  illustrates  the  application  of  modal  coordinate  methods  to  some  representative 
problems  in  Civil  Engineering.  The  three  problems  considered  illustrate  a  range  of  solution 
behaviors,  as  well  as  different  amounts  of  success  that  can  be  encountered  when  a  reduced 
coordinate  analysis  is  used.  The  first  problem  involves  a  dynamic  (blast)  problem,  while  the 
second  is  a  time-dependent  model  that  is  diffusive  in  character.  The  final  problem  is  an  example  of 
how  a  "backward"  reduced  coordinate  scheme  can  encounter  difficulties  in  certain  settings. 

Boundinq-Surface  Soil  Dynamics 

The  first  application  considered  involves  the  dynamic  response  of  a  building  founded  in  relatively 
soft  soil.  The  building  is  loaded  by  a  blast  for  a  duration  of  two  seconds.  The  blast  pressure  starts 
at  zero,  and  grows  linearly  for  those  two  seconds,  reaching  a  maximum  of  four  psi.  For  all  time 
after  the  initial  two- second  blast,  the  pressure  is  taken  as  zero.  Thus,  the  applied  load  acts  for  only, 
the  initial  two  seconds,  and  the  soil-structure  system  responds  in  free-vibration  afterward. 

The  geometry  of  the  problem  is  shown  in  Figure  3,  and  the  appropriate  boundary  conditions  are 
illustrated  in  Figure  4  (this  figure  also  shows  the  location  of  two  particular  elements  that  are  used 
to  monitor  the  behavior  of  the  analysis).  The  building  is  modelled  as  an  isotropic,  linear-elastic 
continuum  in  plane  strain,  which  is  appropriate  for  a  "shear-wall"  building  (where  infilling  panels 
contribute  significantly  to  the  strength  of  the  building’s  lateral-force  resisting  system).  The 
material  properties  for  the  building  are  given  by: 


Ebidg  =  1.0  X  IO7  psf  Vbidg  =  0.18  Tfbidg  =  75  pcf 


Figure  3:  Geometry  for  the  Soil  Dynamics  Problem 
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Figure  4;  Boundary  Conditions  for  the  Soil  Dynamics  Problem 

The  soil  material  is  characterized  using  bounding-surface  plasticity  for  clays  [4,5, 6,7],  and  the 
material  parameters  for  the  soil  are  given  in  Figure  5.  The  clay  soil  is  idealized  as  slightly 
overconsolidated,  with  an  excess  pressure  equal  to  1250  psf  (i.e.,  as  if  the  soil  deposit  had  been 
consolidated  sometime  in  the  past  under  an  additional  10  foot  overburden  with  the  same  unit 
weight). 


Vsoil  —  0.2 

Ysoii  =  125  pcf 

X  =  0.075 

K  =  0.010 

Pi  =  6.35  psf 

Patm  =  2117  psf 

Mc=  1.350 

Me^Mc  =  2/3 

Rc=  3.050 

Re/Rc  =  0.560 

Ac=  0.175 

Ae/Ac  =  0.850 

T=  0.010 

C=  0.485 

S=  1.000 

m=  0.020 

H2=  1.900 

Hc=  2.000 

He/Hc=  0.875 

Figure  5:  Clav  Material  Parameters 

The  initial  state  of  the  clay  under  the  building’s  load  was  obtained  by  solving  the  quasistatic 
problem  of  the  soil- structure  system  loaded  by  its  own  weight.  The  resulting  soil  pressures  were 
used  as  the  initial  stress  slate  for  the  elements  in  the  clay  layer. 

The  analysis  of  the  building  was  carried  out  using  the  Newmark-Newton  method  of  [10].  The 
temporal  discretization  involved  100  time  steps,  each  of  length  0.1  second.  The  use  of  a  longer 
time  step  resulted  in  more  "sub-stepping"  by  the  temporal  integration  scheme,  which  led  to  greater 
computational  cost.  Both  the  direct  (unreduced)  formulation  and  the  reduced  coordinate  approach 
used  incremental  displacements  as  pnmary  unknowns,  since  the  bounding  surface  plasticity  model 
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involves  an  incremental  stress-strain  relation.  In  the  reduced  formulation,  results  were  obtained 
for  2,  4,  and  8  Lanezos  vectors.  In  each  case,  these  vectors  were  generated  with  the  spatial 
distribution  of  the  blast  load  used  as  the  starting  vector  for  the  Krylov  Sequence.  This  choice 
corresponds  to  choosing  the  initial  mode  to  correspond  to  a  pseudo-static  response  to  the  applied 
blast  force.  Adding  vectors  beyond  a  four-dimensional  subspace  did  not  appreciably  increase  the 
accuracy  of  the  analysis,  but  began  to  increase  the  computational  effort.  In  this  problem,  the 
reduced  algorithm  typically  used  approximately  one-third  as  much  computer  time  as  did  the 
unreduced  problem  (even  though  the  code  used  for  the  reduced  algorithm  is  heavily  instrumented 
in  order  to  develop  and  evaluate  the  projection  solution  scheme,  while  the  program  for  the  direct 
formulation  is  a  much  more  efficient  "production"  code).  Because  of  this  fact,  and  also  because 
this  problem  is  still  relatively  small  in  terms  of  the  size  of  the  equation  set,  the  computational 
advantage  of  the  reduced  algorithm  is  somewhat  conservative.  It  is  expected  that  the  reduced 
scheme  will  give  even  greater  efficiencies  on  larger  problems,  especially  those  resulting  from  three- 
dimensional  models. 

The  shapes  of  the  displacement  fields  for  the  first  four  Lanezos  vectors  ("modes")  are  shown  in 
Figure  6.  These  patterns  of  displacement  are  magnified  a  thousand  times  because  the  normalization 
of  these  modes  with  respect  to  the  mass  matrix  yields  actual  displacement  components  that  are  on 
the  order  of  thousandths  of  a  foot.  As  may  be  noticed  in  Figure  6,  the  shape  of  the  displaced 
building  in  the  first  mode  appears  to  incorporate  some  rigid  rotation  of  the  building  under  the 
applied  load,  and  the  next  two  modes  clearly  demonstrate  bending  behavior  of  the  building. 


Magnification  =  lOOOx 


Figure  6:  Lanezos  Modes  for  the  Soil  Dynamics  Problem 
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Some  representative  results  of  the  analyses  are  shown  in  Figures  7  through  11.  In  Figures  7  and 
8,  the  horizontal  displacements  at  the  bottom  and  top  of  the  building  are  plotted  as  functions  of 
time.  It  is  clear  from  these  graphs  that  the  reduced  coordinate  analyses  are  producing  reasonable 
approximations  to  the  overall  behavior  of  the  mechanical  system,  but  that  the  direct  scheme  predicts 
the  permanent  set  of  the  building  (and  the  associated  residual  stresses  in  the  clay)  much  more 
accurately.  This  is  because  the  modal  analyses  involve  projection  vectors  that  do  not  include  large 
plastic  effects,  and  thus  cannot  accurately  represent  the  permanent  strains  due  to  the  blast  loading. 
In  addition,  the  reduced  analyses  underestimate  the  amount  of  energy  dissipation  of  the  nonlinear 
material  model.  In  fact,  if  the  nonlinear  convergence  tolerance  is  relaxed  for  the  reduced  algorithm, 
the  response  approaches  that  of  a  linear  problem,  with  little  or  no  losses  due  to  plasticity. 

It  is  apparent  from  these  results  that  problems  governed  by  inelastic  effects  will  require  special 
attention  from  any  successful  modal  analysis  scheme.  In  particular,  some  means  of  either  updating 
the  modes  to  include  plastic  deformation,  or  of  increasing  the  size  of  the  projection  subspace  will 
be  required  for  high  accuracy  from  this  family  of  approximate  solution  techniques.  Given  the  very 
slight  increase  in  accuracy  observed  in  this  problem  when  the  number  of  modes  is  increased,  it 
appears  that  modal  updating  schemes  are  the  most  appropriate  path  for  future  work  in  this  area. 


Figure  7:  Horizontal  Displacements  at  the  Bottom  of  the  Building 
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Figures  9,  10,  and  11  depict  representative  building  and  soil  stresses  for  the  elements  highlighted 
in  Figure  4.  These  stresses  show  better  agreement  between  direct  and  reduced  results  than  the 
displacements  displayed  in  Figures  7  and  8.  As  in  the  comparison  of  displacements,  the  addition 
of  more  modes  did  not  produce  significantly  better  results,  which  reinforces  the  view  that  some 
means  of  updating  a  small  number  of  modes  is  probably  a  more  efficient  reduced  coordinate 
strategy  than  increasing  the  size  of  the  projection  subspace. 


This  example  demonstrates  that  the  reduced  scheme  can  be  used  to  obtain  reasonable 
approximations  to  difficult  problems  in  dynamics,  and  that  such  reduced  coordinate  algorithms 
warrant  further  study.  Even  with  the  present  underestimation  of  the  plastic  effects  in  this  problem, 
the  reduced  approach  is  useful  as  a  less  expensive  tool  for  a  preliminary  analysis.  In  addition,  its 
memory  requirements  are  extremely  small  compared  to  the  direct  method,  since  the  global 
stiffness,  tangent  stiffness,  and  mass  matrices  do  not  need  to  be  formed,  stored,  or  factored.  The 
next  example  problem  demonstrates  a  nonlinear  time-dependent  problem  where  these  memory 
needs  are  very  important. 
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Two  Lanczos  Vectors 
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Figure  1 1 :  Soil  Shear  Stresses  at  Base  of  Building 


Another  important  time-dependent  phenomena  involving  soils  involves  the  flow  of  fluids  through 
soil  voids.  Examples  of  this  problem  are  soil  consolidation  (where  the  fluid  is  porewater),  gas 
migration  (where  the  fluid  is  air  or  methane),  and  contaminant  transport  (where  a  contaminating 
fluid  is  carried  through  the  soil  voids  by  the  motion  of  groundwater).  All  of  these  are  important 
"real-world"  problems  of  considerable  interest  to  engineers,  and  all  can  be  approached  from  a 
reduced  coordinate  standpoint.  Methods  for  practical  solution  of  soil  consolidation  problems  have 
been  addressed  in  [7],  and  contaminant  transport  problems  are  governed  by  a  non-self-adjoint 
boundary- value  problem  that  is  an  object  of  present  research  by  the  authors  (recall  that  much  of  the 
theoretical  framework  presented  in  the  last  chapter  involved  self-adjoint  problems).  The  problem 
of  gas  migration  through  soils  is  governed  by  a  generalization  of  the  familiar  diffusion  equation, 
and  is  a  good  candidate  for  a  reduced  coordinate  formulation. 
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In  this  case,  the  fluid  is  methane  gas  produced  by  the  decomposition  of  solid  waste  in  a  municipal 
landfill,  and  the  porous  medium  includes  both  soil  and  the  waste  itself,  but  the  basic  principles  and 
solution  process  for  this  problem  are  readily  applicable  to  many  other  problems  involving  diffusive 
phenomena.  In  the  setting  of  gas  production,  this  model  can  be  taken  as  a  "zeroth"  order 
approximation  to  a  very  complicated  problem  that  is  of  urgent  importance  in  many  urban  areas. 
Many  of  the  simplifications  made  are  somewhat  crude,  but  they  serve  as  a  good  starting  point  for  a 
more  refined  physical  and  mathematical  analysis.  The  primary  interest  in  this  problem  from  the 
standpoint  of  this  reduced  coordinate  research  is  that  the  gas  model  produces  an  interesting 
nonlinear  problem,  and  the  reduced  scheme  promises  to  allow  reasonable  approximate  solutions  to 
be  found  using  microcomputers  such  as  those  that  can  be  found  in  many  small  engineering  firms. 


The  governing  equation  for  two-dimensional  flow  in  a  porous  medium  is  given  by  the  diffusion 
equation; 


c 


a 

57 


+  bu  =  f 


The  functions  c,  kx  ,  ky ,  and  b  are  material  parameters,  u(x,y)  is  the  solution  (pressure,  in  this 
case),  and  f  is  a  source  term,  which  can  include  singular  behavior  such  as  point  sources.  The 
coefficient  c  represents  the  "storage"  capacity  of  the  medium,  the  functions  kx  and  ky  represent  the 
medium's  ability  to  conduct  fluid,  and  the  term  b  reflects  the  convection  of  fluid  in  the  direction 
perpendicular  to  the  plane  (i.e.,  leakage  to  the  atmosphere).  This  equation  is  a  first-order 
differential  equation  in  time,  and  similar  equations  govern  a  wide  variety  of  diffusive  problems. 

In  the  cciuext  of  the  gas  production  problem,  the  source  function  f  depends  on  time,  since  the 
bacterial  decomposition  process  that  produces  the  gas  is  dependent  on  temperature  and  water 
supply,  both  of  which  may  vary  with  the  season.  In  addition,  because  there  is  a  period  of  time 
between  the  placement  of  the  landfill  and  the  onset  and  development  of  microbial  decomposition, 
the  time  dependence  of  the  source  term  is  modelled  by  separating  this  term  into  two  parts: 

f(x,y,t)  =  h(t)  s(x,y) 

The  function  h(t)  represents  a  "history"  function  whose  long-range  average  is  unity,  and  this 
function  scales  the  spatial  source  term  s(x,y).  The  history  function  h(t)  used  in  this  example 
problem  is  shown  in  Figure  12. 
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Figure  12:  History  Function  for  Gas  Production 


The  domain  of  the  problem  is  actually  three-dimensional,  but  a  horizontal  plane  in  the  landfill  is 
modelled  as  a  simple  two-dimensional  approximation.  It  is  assumed  that  the  soil  cover  above  the 
landfill  allows  some  migration  of  the  gas  in  a  direction  perpendicular  to  the  plane,  and  that  this 
migration  is  dependent  on  the  pressure  in  the  landfill,  since  large  enough  pressures  should  allow 
the  gas  to  breach  the  soil  cover.  In  this  problem,  this  variation  in  perpendicular  permeability  of  the 
cover  is  modelled  using  a  nonlinear  convection  function  b  =  b(x,  y,  u)  defined  by: 

b(x,y,u)  =  bo(x,y)  p(x,y)  u2 

This  quadratic  convection  coefficient  represents  a  soil  cover  that  "heals"  after  being  subjected  to 
large  pressures.  Many  other  choices  are  possible,  and  the  most  appropriate  representation  for  this 
term  is  still  under  study.  The  material  parameters  used  for  the  model  are  tabulated  in  Figure  13, 
and  the  two  dimensional  domain  is  shown  in  Figure  14  . 
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Mv  +  Kp  =  f 
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where  v  =  9p/9t  is  the  time  rate  of  change  of  the  vector  p  of  nodal  pressures.  Formulae  for  the 
mass  M  and  stiffness  K  can  be  found  in  any  standard  Finite  Element  text,  though  it  should  be 
noted  that  the  stiffness  K  depends  on  p,  so  this  problem  is  governed  by  a  nonlinear  system. 

The  required  temporal  integration  scheme  can  be  developed  by  expressing  the  pressure  p  at  time 
tn+i  in  terms  of  the  "velocity"  v  at  tn+i  and  terms  evaluated  at  time  tn: 


pn+l  =  pn  +  (1  -  a)vnh  +  ttvn+lh 


(h  =  tn+l  -  tn) 


This  relation  for  pn+l  in  terms  of  vn+i  can  be  inverted  to  cast  the  pressure  p  as  the.  primary 
independent  variable; 

1  ,  ,  (1-a) 

Vn+l  =  —  (Pn+l  -  Pn  ) - 

ah  a 

A  "predictor"  step  can  be  formed  by  using  this  result  to  write  the  governing  matrix  equation  at  the 
end  of  the  time  step  (t  =  tn+i)  entirely  in  terms  of  the  pressure  vector: 

Keffpn+l  =  feff 

where  =  K  +  —  M 
ah 

r  r  w  r  1  (1-a)  1 

feff  =  fn+1  +  M  [  — Pn+ - V„  ] 

ah  a 

A  Newton  "corrector"  scheme  can  be  derived  in  this  case  (where  M  and  f  are  not  functions  of  the 
pressure  p)  by  defining  a  residual  vector  r: 

r(p)  =  K(p)  p  +  Mv  -  f(p)p 


The  (vector)  zeros  of  this  residual  can  be  found  by  solving  the  relation: 


r'((i-i)p)  [(i)P  -  (i-i)p]  =  -r((i.i)p) 

...  r  9K  3f 

With  r(p)  =  [-^p  +  K]  +  M^-^ 

r  dK  1  . .  9f 

=  [-5— p  +  K]+  —  M-^ 
^  ah  ^ 
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(where  i  is  the  iteration  counter)  until  convergence  to  a  fixed  point  is  obtained.  In  the  gas- 
production  problem,  the  quadratic  nonlinearity  in  b(x,y)  leads  to  a  symmetric  tangent  stiffness 
[(9K/3p)p  +  K].  More  complicated  nonlinearities  (e.g.,  those  affecting  the  conductivity  terms) 
may  lead  to  a  non-symmetric  tangent  stiffness. 

This  "predictor-corrector"  algorithm  was  used  to  solve  the  gas-production  problem  for  three  cases: 

(1)  A  direct  formulation  involving  the  mesh  of  Figure  14 

(2)  A  direct  formulation  for  a  refined  mesh  obtained  by  subdividing  each  element  from  analysis  (1) 

into  four  elements. 

(3)  Reduced  solutions  involving  1  to  4  Lanezos  vectors  and  the  coarse  mesh  of  analysis  (1). 

The  four  Lanezos  vectors  ("modes")  used  in  analysis  (3)  are  shown  in  Figure  15. 


Figure  15:  Lanezos  Modes  for  Gas  Flow  Problem 

Results  from  the  various  analyses  are  shown  in  Figures  16,  17,  and  18.  Plots  of  pressure  against 
time  are  found  in  Figures  16  and  17  for  the  two  sampling  points  E  and  F  shown  on  Figure  13. 
The  "exact"  solution  (black  line)  is  the  curve  for  the  coarse  mesh  —  results  for  this  mesh  are 
practically  identical  to  those  from  the  refined  mesh,  leading  to  the  conclusion  that  this  solution  has 
essentially  converged.  The  modal  analyses  that  are  not  shown  coincide  with  this  "exact"  solution 
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and  are  not  plotted.  Figure  18  shows  the  distribution  of  pressure  at  the  end  of  the  analyses.  The 
modal  results  are  excellent  approximations  to  the  "exact"  solution  even  after  60  time  steps  (3 
years). 


0.0 


3.0 
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Figure  18;  Final  Solution  Contours  for  Gas  Flow  Problem 


This  example  is  clearly  a  successful  application  of  a  reduced  coordinate  method  to  an  important 
problem  in  Civil  Engineering.  The  solution  efficiencies  encountered  here  promise  that  three- 
dimensional  problems  of  this  type  should  be  solvable  even  on  low-cost  microcomputers.  The 
actual  behavior  of  this  solution  depends  almost  entirely  on  the  "low-frequency"  modes  (recall  that 
there  are  no  dynamic  effects  here),  and  the  original  mode  shapes  capture  the  solution  behavior  even 
after  substantial  nonlinearities  have  developed.  This  is  a  typical  result  for  a  problem  for  a  self 
-adjoint  (in  the  linear  case)  problem  where  the  data  primarily  excites  the  lowest  frequencies.  The 
next  example,  however,  demonstrates  that  this  excellent  agreement  between  reduced  coordinate 
approximation  and  direct  formulation  is  not  always  observed. 
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Soil  Dynamics  Incorporating  Porewater  Effects 

The  response  of  soils  of  low  permeability  (such  as  clays  and  fine  silts)  to  dynamic  loads  is  largely 
governed  by  the  rate  at  which  porewater  trapped  in  the  soil  voids  can  flow.  The  soil-water  system 
thus  acts  like  a  large  sponge  in  that  the  deformations  of  the  solid  phase  are  influenced  by  the  rate  at 
which  the  fluid  phase  responds  to  changes  in  pressure.  Thus  any  completely  realistic  model  of 
structural  behavior  for  cohesive  soils  will  have  to  consider  the  effects  of  pore  pressure.  One  of  the 
most  important  of  these  effects  is  that  a  relatively  impermeable  soil  behaves  as  a  nearly 
incompressible  medium,  which  often  necessitates  using  low-order  (reduced)  integration  on  various 
volumetric  energy  terms  in  order  to  prevent  the  Finite  Element  solution  from  "locking-up"  (i.e., 
seriously  underestimating  the  magnitude  of  displacements). 

For  the  problem  of  soil  consolidation,  the  structural  response  of  the  soil  mass  is  quasistatic,  and 
inertial  effects  can  be  neglected.  This  class  of  time-dependent  problems  is  developed  in  [7],  where 
a  Finite  Element  analysis  for  consolidation  of  clay  soils  is  presented.  The  soil  dynamics  model 
used  in  this  third  and  last  example  uses  the  basic  equations  of  this  Finite  Element  soil  code,  with 
inertial  terms  added  to  account  for  dynamic  effects.  The  behavior  of  porewater  is  modelled  by 
inclusion  of  additional  displacements  at  each  node,  which  represent  the  movement  of  the  porewater 
relative  to  the  soil  skeleton.  A  detailed  derivation  of  this  entirely  displacement-based  formulation 
of  the  soil  dynamics  problem  can  be  found  in  [17].  In  this  reference,  it  is  mentioned  that  the 
formulation  may  suffer  from  "hourglassing",  which  is  an  oscillatory  phenomenon  that  erodes  the 
quality  of  the  solution.  This  hourglassing  response  arises  from  the  introduction  of  spurious  low- 
energy  deformation  modes  produced  as  an  unwanted  by-product  of  the  low-order  integration 
terms.  Technically,  these  hourglassing  effects  are  termed  spurious  singular  modes,  because  they 
are  unwanted  (and  incorrect)  solution  "shapes"  that  correspond  to  near- zero  eigenvalues  of  the 
Finite  Element  equation  set. 

The  effect  of  these  singular  modes  is  shown  in  Figure  19,  which  is  a  plot  of  two  Lanezos  vectors 
obtained  from  a  backward  modal  projection  scheme  (compare  with  Figure  6).  The  problem  that  is 
modelled  is  the  same  as  the  first  example  problem  of  this  chapter,  but  generalized  in  the  manner  of 
the  last  paragraph  to  include  porewater  flow.  The  direct  (unreduced)  solution  gives  reasonable 
results  (which  are  similar  to  those  of  Figures  7  through  1 1),  but  the  results  of  the  reduced 
coordinate  formulation  for  this  case  are  completely  incorrect,  and  underestimate  the  solution 
displacements  by  several  orders  of  magnitude.  The  two  Lanezos  vectors  shown  give  the  reason 
for  this  complete  degradation  of  solution  quality;  the  dynamic  response  of  the  soil-building 


Modal  Projcciion  Methods 


Page  24 


Structure  does  not  even  appear  in  these  modes  used  for  the  reduced  solution,  so  that  the  reduced 
coordinates  cannot  account  for  any  movement  of  the  structural  system.  The  problem  is  that  the 
spurious  hourglassing  shapes  appear  in  the  problem  as  extremely  low-frequency  vibrational 
modes,  which  are  then  detected  and  used  by  the  backward  projection  scheme  to  the  exclusion  of 
the  correct  (i.e.,  non-spurious)  modes.  The  ability  of  backward  projection  methods  to  identify  the 
low-frequency  response  of  a  structural  system  has  been  obviated  by  the  fact  that  the  lowest 
frequencies  of  this  Finite  Element  model  are  entirely  spurious! 


i 


Figure  19:  Sample  Hourglass  Modes  for  Soil  Dynamics  Problem 

This  result  could  have  been  predicted  by  the  analysis  of  Chapter  2,  where  it  was  shown  that  a 
structural  system  possessing  near-instabilities  would  cause  difficulties  for  a  backward  projection 
method.  This  problem  (in  its  present  form)  would  be  better  solved  by  a  forward  projection  scheme 
such  as  a  preconditioned  forward  Lanczos  method,  which  would  converge  from  the  high- 
frequency  end  of  the  spectrum  and  thus  avoid  the  spurious  hourglassing  modes. 
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A  better  approach  would  involve  n2l  using  models  with  pathological  behavior  such  as  these 
singular  modes.  The  problem  with  hourglassing  can  be  removed  by  using  a  mathematically  more 
correct  mixed  formulation  and  a  better  choice  of  element  interpolants,  or  by  adding  fictitious 
viscosity  to  the  model  for  the  purpose  of  "hourglass  control".  This  example  is  intended  to  show 
that  "backward"  reduced  coordinate  methods  cannot  be  blindly  applied  to  arbitrary  types  of  Finite 
Element  problems,  but  only  to  those  cases  where  the  solution  is  primarily  dependent  on  a  relatively 
small  number  of  lower-frequency  modes  of  response. 

These  examples  demonstrate  the  notion  that  reduced  coordinate  models  are  a  useful  tool  for  cutting 
many  large  intractable  problems  down  to  a  more  manageable  size,  but  this  approach  should  be  used 
with  caution  where  pathological  problems  with  the  formulation  are  known  to  exist.  Modal  analysis 
procedures  are  thus  one  more  computational  tool  for  solving  large  problems,  but  should  not  be 
viewed  as  a  panacea  for  curing  all  the  difficult  problems  in  computational  mechanics. 
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Extensions 

Suggestions  for  Further  Study 

The  reduced  coordinate  scheme  introduced  in  [10]  and  applied  in  this  document  can  be  seen  to  be  a 
useful  solution  tool  for  many  types  of  engineering  problems,  but  there  is  clearly  more  work  that  is 
needed  to  make  this  technique  more  reliable,  robust,  and  efficient.  A  few  of  the  areas  where  future 
development  and/or  application  is  warranted  include: 

(1)  Application  of  the  method  to  framed  structures.  These  systems  are  typically  very  dependent 
upon  low-fi^uency  response,  and  are  thus  likely  candidates  for  efficient  application  of  the 
proposed  reduced  coordinate  scheme.  Present  research  by  the  authors  involves  application  of 
these  modal  methods  to  two-  and  three-dimensional  frames  which  incorporate  nonlinear 
material  effects  through  the  use  of  bounding-surface  plasticity  models  relating  moment  and 
curvature.  Preliminary  efforts  in  this  area  have  demonstrated  that  the  reduced  coordinate  model 
is  an  appropriate  tool  to  reduce  the  size  of  problems  of  this  class. 

(2)  Application  of  the  method  to  soil  dynamics  problems  involving  porewater  effects.  Although 
the  naive  implementation  of  the  reduced  coordinate  model  for  the  last  problem  of  the  preceding 
chapter  showed  poor  results,  it  appears  that  this  method  can  be  useful  in  this  case,  as  long  as  it 
is  applied  to  a  formulation  that  does  not  exhibit  spurious  singular  modes.  Attention  is  being 
focused  on  mixed  Finite  Element  models  that  satisfy  consistency  conditions  on  displacement 
and  pressure  interpolants  designed  to  inhibit  the  existence  of  these  spurious  modes.  In 
addition,  generalization  of  the  theoretical  framework  of  the  projection  scheme  to  accommodate 
non-self-adjoint  problems  (such  as  those  arising  from  bounding-surface  characterizations  for 
sands)  is  an  object  of  continuing  research. 

(3)  Development  of  schemes  for  updating  modes  to  account  for  inelastic  effects.  This  research  is 
proceeding  for  both  continuum  and  framed  structures,  and  is  concentrated  on  efficient  modal 
update  schemes  that  do  not  result  in  any  truncation  errors  from  the  change  in  projection  basis. 
Since  the  Finite  Element  method  itself  is  a  projection  method,  one  promising  avenue  of 
research  is  to  view  the  modal  update  procedure  as  an  analog  of  a  Finite  Element  mesh  rezoning 
problem  (where  many  proposed  schemes  have  been  developed).  Mapping  the  solution  from 
one  modal  basis  to  another  without  truncation  error  may  be  affected  in  a  manner  similar  to  a 
mesh  remapping  algorithm.  This  topic  is  the  subject  of  continuing  research  by  the  authors. 
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The  reduced  coordinate  schemes  applied  in  this  report  show  considerable  potential  to  solve  a 
variety  of  important  engineering  problems.  Placed  in  the  context  of  projection  methods,  they  can 
be  seen  as  part  of  a  unified  array  of  methods  that  includes  the  Conjugate  Gradient  Method,  the 
Lanczos  Algorithm,  and  modal  superposition  schemes.  All  of  these  approaches  share  the  common 
characteristic  that  each  may  be  used  to  make  certain  large  intractable  problems  solvable  on  a  variety 
of  computers,  and  thus  extend  the  range  of  problems  that  can  be  solved.  The  examples  presented 
demonstrated  a  spectrum  of  success  in  the  application  of  the  backward  modal  analysis  schemes 
ranging  from  excellent  agreement  to  near-total  failure.  On  non-pathological  problems  typical  of 
those  encountered  in  practice,  the  reduced  scheme  performed  well  in  that  essential  phenomena  were 
captured  while  decreases  in  computer  expense  were  observed.  On  the  "ill-behaved"  problem 
presented,  the  modal  scheme  gave  poor  results,  but  the  source  of  the  difficulty  was  seen  to  be  a 
shortcoming  in  the  formulation  of  the  unreduced  model,  and  the  poor  results  could  be  alleviated  by 
a  more  appropriate  choice  of  projection  scheme. 

In  summary,  the  reduced  coordinate  method  has  demonstrated  its  utility  on  a  number  of  important 
problem  in  Civil  Engineering,  and  warrants  further  work  to  extend  the  range  of  its  usefulness. 
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